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Abstract. Spin tunneling in molecular magnets controlled by dipole-dipole interactions (DDI) in the dis- 
ordered state has been considered numerically on the basis of the microscopic model using the quantum 
mean-field approximation. In the actual case of a strong DDI, coherence of spin tunneling is completely lost 
and there is a slow relaxation of magnetization, described by t^^^ at short times. Fast precessing nuclear 
spins, included in the model microscopically, only moderately speed up the relaxation. 



PACS. 75.45.-|-j Macroscopic quantum phenomena in magnetic systems - 
resonances and relaxations - 75.75.Jn Dynamics of magnetic nanoparticles 



76.20.-|-q General theory of 



1 Introduction 

Molecular magnets (MM) built of molecules with a large 
effective spin, such as S* = 10 in Mni2 and Fes, attract 
attention by their bistability resulting from a large uni- 
axial anisotropy D |1|2] that creates the energy barrier 
DS^ ~ 67K for spin rotation (for a review, see |3l4l5j ). 
The most spectacular finding on molecular magnets is 
resonance spin tunneling |6l7l8j that occurs when spin 
energy levels on the different sides of the barrier match 
each other. Another striking phenomenon is propagating 
fronts of magnetic burning or deflagration |9I10) . a self- 
supporting process driven by the energy release that con- 
trols thermally activated overbarrier relaxation As 
further development, fronts of spin tunneling controlled by 
the dipolar field |12I13] and a combined quantum-thermal 
theory of magnetic deflagration [Tl] have been proposed. 

Magnetic molecules in MMs form a crystal lattice (bo- 
dy-centered tetragonal for Mni2 Ac) . As magnetic cores of 
the molecules are shielded by organic ligands, there is no 
exchange interaction between the molecules in the crystal, 
thus the dipole-dipole interaction (DDI) is dominating. 
There is an evidence of dipolar ordering below 1 K in Mni2 
Ac and other MMs [15116117] . dynamics of which is due 
to spin tunneling. Relaxation of spin states in molecular 
magnets due to interaction with the environment should 
be a collective process, as many molecular spins are inter- 
acting with the same phonon or photon modes. Theories 
of photon and phonon superradiance in MMs have been 
proposed in Refs. |18ll9j . 

In the limit of low temperatures, processes of thermal 
activation die out and the only way for the molecular spin 
to cross the barrier is spin tunneling. Coherence of spin 
tunneling between the two degenerate ground states in 
molecular magnets is difficult to observe because the en- 
ergy bias due to DDI is much greater than the tunnel 
splitting of spin states A. As a result, most of molec- 
ular spins are strongly biased and cannot tunnel. Fur- 
thermore, tunneling of a spin changes the dipolar field 
on other spins, prohibiting or allowing them to tunnel. 
Thus low-temperature spin relaxation in molecular mag- 
nets is a complicated collective process that results in a 



non-exponential slow relaxation. The latter was initially 
observed as magnetization relaxation from the saturated 
state in Fes EQ] and Mni2Ac [21] that followed the Vi 
law. Theoretically this behavior had been explained by 
spreading of the initially uniform dipolar field in samples 
of ellipsoidal shape in the course of relaxation [22 23i24] . 
This spreading gradually drives more and more spins off 
resonance and slows down the relaxation. The situation is 
much less clear in the realistic case of MM crystals of ar- 
bitrary shape where the dipolar field is non-uniform from 
the beginning. 

Tunneling of most of the spins being hampered by a 
large dipolar bias leads to the quest for a mechanism that 
could accelerate spin transitions. It was suggested that 
nuclear spins provide a fast fluctuating bias on magnetic 
molecules that is much larger than A and brings a larger 
number of molecules on resonance, allowing them to tun- 
nel [35]. However, combining formulas of Ref. [32] leads 
to the final result for the relaxation rate that does not 
depend of nuclear spins [24] . The controversy penetrated 
experimental literature as well, where the formula for the 
^/i relaxation rate (see, e.g., Eq. (2) of Ref. [35] ) 



A/2 ^ A^P{^H)/tl, 



(1) 



is used, P{£,h) being the distribution of the energy bias 
on magnetic molecules. As the latter is mainly due to the 
DDI, there is no place for nuclear spins in this formula. 
With no external bias one has P{£.h) ^ ^/Ed and [33] 



A/2 - A^ihEL 



(2) 



where Ed is the dipolar energy {Eo/ks is 0.067 K for 
Mni2 Ac and 0.126 K for Feg). With A/ke ^ 10"^ K for 
Fes at zero field, the factor suppressing spin tunneling is 
very large: Ed/ A '-^ 10^. 

Further experimens showed that the y/t relaxation also 
follows a small abrupt change of the bias field in the disor- 
dered state of MMs [25126) . Here it is difficult to propose 
an analytical approach, and numerical methods had been 
applied. Refs. |27I28| state that the ^/t relaxation law is 
not universal, and for the body- and face-centered lattices 
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the relaxation exponent is 0.7 rather than 1/2. These re- 
sults were criticized in Refs. |29I30I31] as finite-size effect 
and effect of fitting the data over a too long time inter- 
val, and \/t relaxation law was obtained for all lattices 
at times short enough. However, Refs. [32133) refute the 
criticism, insisting on the non-universal behavior with the 
relaxation exponent of up to 0.73 for the face-centered lat- 
tice. At larger times, the magnetization relaxation follows 
a stretched exponential [55] . 

Recently dipolar-controlled relaxation together with 
dipolar ordering has been observed in dynamic suscep- 
tibility experiments on an Er-based MM |34) . 

Most of published numerical work use Monte Carlo 
simulations based on the "tunneling window" concept, ac- 
cording to Ref . [33] . Spins are being checked one after the 
other, and if the bias on a spin is within the preset tunnel- 
ing window, the spin is allowed to fiip. The ensuing change 
of the dipolar fields is taken into account when checking 
next spins. The tunneling window in these simulations had 
been set to the amplitude of the fiuctuating bias due to 
nuclear spins that is much larger than A. Although this 
approach captures the essential physics of the DDI and 
leads to the results qualitatively similar to what is seen 
in experiments, it is oversimplified and postulates the role 
of nuclear spins instead of describing their action dynam- 
ically. The ensuing relaxation rates are roughly propor- 
tional to the tunneling window due to nuclear spins, that 
contradicts Eq. On the other hand, the absolute value 
of the relaxation rate cannot be found by Monte Carlo sim- 
ulations because Monte Carlo steps cannot be quantified 
in terms of real time without additional work. The recent 
work [35[ studying low-temperature relaxation in molecu- 
lar magnets out of saturation with Monte Carlo and then 
rate equations, takes nuclear spins into account in a differ- 
ent, although also phenomenological way. Ref. |36j consid- 
ers a one-molecule model with a stochastic random-walk- 
type bias, also obtaining the \/t law. Ref. [371 considers 
the influence of nuclear spins on tunneling in the absence 
of DDI. 

Theoretical methods mentioned above replace the orig- 
inal quantum-mechanical model by a more tractable sto- 
chastic model. Clearly, the ultimate solution of the prob- 
lem should use the many-body Schrodinger equation. This 
approach has been developed for a "cental" system con- 
sisting of one spin or several interacting spins, coupled to 
a sea of other spins, e.g., nuclear spins. The most efficient 
numerical method is based on the Chebyshev expansion 
[38139] . Since the number of coefficients in the expansion 
of the wave function of the system over a basis increases 
exponentially with the system size, calculations are fea- 
sible for a number of spins up to 20 and require a su- 
percomputer. Decoherence and approach to equilibrium 
could be established with this method for small systems 
[40141142] . However, solution of the full Schrodinger equa- 
tion for the system under consideration is problematic be- 
cause the dipole-dipole interaction is long-ranged and the 
size of the system has to be much larger. 

The aim of the present work is to provide a microscopic 
description of the dipolar-controlled tunneling relaxation 



in molecular magnets on the basis of the real dynamics 
of molecular and nuclear spins, and to clarify the long- 
standing controversy of the role of the latter. The method 
consists in numerically solving the system of equations of 
motion for molecular spins at equilibrium in the quan- 
tum mean-field (Hartree) approximation and measuring 
the time-dependent auto-correlation function. The latter 
with the help of the linear-response theory and other re- 
lations can be proven to be proportional to the magneti- 
zation relaxation function measured in experiments. This 
approach is presented in Sec. [51 The numerical procedure 
is described in Sec. [31 while Sec. HI contains numerical re- 
sults for the Mni2 Ac and simple cubic lattices in the 
absence of an external bias. In Sec. [S] a stronger-than- 
expected influence of an external bias on spin relaxation 
is shown. Section [7| introduces nuclear spins microscopi- 
cally via time-dependent bias fields acting on the molec- 
ular spins. The results show a moderate speed-up of the 
relaxation due to nuclear spins in the case of their fast 
precession. In the concluding section, the results obtained 
and their possible generalizations are discussed. 

2 The model and equations of motion 

At low temperatures, only the lowest doublet of the Hamil- 
tonian of a magnetic molecule is populated, so that one 
can use the two-state model described by the pseudospin 
1/2 

i i ij 

Here cr is the Pauli matrix, A is tunnel splitting, Wq = 
2Sgfj,BBz is the external bias and Vij is the dipole-dipole 
interaction 

^ , , 3 cos 9ij — 1 
Vij^Eocpij, (Pij=vo 3 , (4) 

where 

Ed = {SgfiB)Vvo (5) 

is the dipolar energy {Eo/kB — 0.067K for Mni2 Ac), 
Wq is the unit-call volume, is the distance between the 
molecules at sites i and j , and 9ij is the angle between 
and the easy axis z. Eq. ([3]) is a kind of transverse Ising 
model with a long-range interaction. 

The Heisenberg equation of motion for spin operators 
<Ti has the form 

a^{t)^a,{t)x fl,, (6) 

where 

(Zie, + W,e,) , = VKo + 2 ^ V.^a,, . (7) 

j 

One can see that Eq. (|6]) describes precession of the spin cr^ 
in an effective field. However, going over from the Heisen- 
berg operators to their quantum-mechanical averages that 
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we are interested in is notrivial because, in general, differ- 
ent spin operators are entangled and their quantum aver- 
ages do not factorize, {<Jia<Jjp) ^ {oia) {(Jjp)- In the ab- 
sence of factorization, one has to resort to the full Schrodin- 
ger equation for a many-spin system that can be solved 
only for a small number of spins |38I39I40I41I42| . To make 
the problem tractable numerically for a larger number of 
spins, one can apply the quantum mean-field or Hartree 
approximation for spins on different sites, 

(o-icO-j/j) =^ {(Tia) {(Tjfi) , j. (8) 

This approximation means that each spin is exposed to the 
quantum-mechanically averaged field from its neighbors, 
that should be justified when the number of interacting 
neighbors is big. For the long-range dipole-dipole inter- 
action, the latter should be the case. After decoupling of 
quantum correlations, the equation of motion for the av- 
erages of spin components has the same form as Eq. (jH]). 
Thus for the ease of notation we replace (<Tj) ^ o-j and 
obtain exactly Eq. now for precessing classical-like 
vectors. 

In the temperature range well above the dipolar order- 
ing temperature ^1 K (T = oo at the scale of the DDI), 
molecular spins are completely disordered and spins of 
different molecules are uncorrelated. This is an additional 
justification for the decoupling in Eq. ([8]). The correlation 
function needed to calculate the spin relaxation function 
thus reduces to the zz component of the autocorrelation 
function 

S^o.p {t) = i {t)a.,p (0) + a.,p (0)a,„ (t)} (9) 
{a, (3 — x,y, z) that has to be averaged over all N spins. 

The autocorrelation function cannot be decoupled because 
this would violate the properties of spin operators. One 
can use the identities 

CzCTz = CTxCx = CTyCTy = 1, (11) 

and 

C^zCTx + CTxCTz = CTzf^y + CTyf^z = (12) 

to obtain the initial conditions for the autocorrelation 
function of Eq. © 

S.apiO)=r' ""^J (13) 

The equation of motion for the autocorrelation function 
can be obtained by differentiating over time, using the 
Heisenberg equation of motion, Eq. ([5]), and decoupling 
i2i since the latter is due to spins on other sites. This 
yields the equation 

^S,,(i) = S,,(i) X f2„ (14) 



where 

S^z{t) - \ {cT,{t)a,M+a,M<^-Xt)) ■ (15) 

Eqs. ^ and (ITil) are coupled and have to be solved to- 
gether. The initial condition for Eq. (jl4l) is 8^2(0) = e^. 
For Eq. dH]), we use the initial condition describing 0-^(0) 
pointing in random directions. Note that the autocorre- 
lation function defined as (Jiz{t)cFiz{Q) is wrong for MMs 
because it describes a system of classical rather than of 
quantum spins. This model is of its own interest and it 
will be briefly addressed in Sec. HI 

With the help of the linear response theory and other 
relations (see Appendix) one can find the spin relaxation 
function measured in the experiment. The latter turns out 
to be just Szz{t) of Eq. (ITU)) . In particular, if a smaU bias 
field Bz is abruptly removed, the ensuing relaxation is 
described by a remarcably simple formula 

IJzit) = (TzMqSzzit), (16) 

where crz,cq — "^Sg^BBz/ {ksT) is the equilibrium spin av- 
erage in the high-temperature limit. Note that only CTz^eq 
contains the temperature, whereas Szz{t) can be calcu- 
lated at T = 00. In other cases, such as a small field Bz 
abruptly applied to an initially unbiased MM at equilib- 
rium, one obtains similar formulas. 

3 Numerical solution 

Numerical solution of Eqs. © and (ITT)) is a serious task 
because there are six dynamic variables per lattice site 
and the DDI is long-ranged. The latter can be dealt with 
by methods based on the fast Fourier transform (see, e.g., 
Ref. [43"). Fortunately, the autocorrelation function self- 
averages over random initial orientations of spins. Thus 
for a very large number of molecular spins the solution 
is practically the same for any realization of the initial 
conditions. However, when the number of spins is not ex- 
tremely large, averaging over several runs with different 
realizations of initial conditions is needed to make the time 
dependence of the autocorrelation function smoother. 

In most of the calculations, the body-centered tetrag- 
onal lattice of Mni2Ac has been used. Some calculations 
have been performed for a simple cubic lattice to check 
universality of the relaxation exponent. Numerically it is 
impossible to deal with the realistic case of very large 
Ed I A because the relaxation becomes too slow and com- 
putation time becomes too long. However, the laige-Eu 
behavior sets in already for Ed/ A > 3, and the results for 
3 > Ed /A > 100 can be scaled to obtain the result for 
any large Ed/ A. 

The sizes of crystals are expressed in units of a for the 
body-centered tetragonal lattice of Mni2 Ac with lattice 
parameters a — b > c. The "external" sublattice con- 
tains NaNbNc spins while the "internal" sublattice con- 
tains {Na - l){Nt - l){Nc - 1) spins, where Na^ La + 1, 
Nb = Lb + I, and Nc = Round(Lc/??) + 1 with ij = c/a. 
The total number of spins is iV = NaNbNc + {Na-l){Nb- 
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l){Nc — !)■ In most calculations, N was above 20000 that 
is comparable with the system sizes in Monte Carlo sim- 
ulations of Ref. [53] and subsequent works. 

The rate of spin tunneling (i.e., precession of the spin 
vector around the x axis) is determined by A. Most of the 
time, however, the bias of the spin is very large, Wi » A, 
and the spin vector is fast precessing around the z axis. 
This fast precession does not contribute to spin tunneling 
and is actually irrelevant. However, it necessitates a very 
small integration step in numerically solving the equa- 
tions of motion, to achieve a sufficient accuracy. For nu- 
merical efficiency, it is absolutely crucial to introduce a 
cutoff of the bias Wi that allows to increase the integra- 
tion step. In these calculations, a cut-off at WVax = 5A 
was introduced by replacing in the equations of motion 
Wi ^ Wi,cff = Wj/Q, a nd A ^ Z\,,eff = A/Q^, where 

Q^ = yiTTWWWF and c = 6. 

Numerical calculations have been performed using Wol- 
fram Mathematica 8.0.4.0 that supports vectorization and 
makes use of Intel's math kernel library (MKL) for num- 
ber crunching. The latter employs threading for standard 
tasks such as solving vectorized systems of ordinary differ- 
ential equations (ODE), that results in a significant speed- 
up on multi-processor machines for a large system size. 
The most efficient algorithm for solving our ODE's proved 
to be Runge-Kutta 4th order with a fixed step size. The 
latter was chosen so that the average deviation of the spin 
length from 1 at the end of the calculation does not exceed 
1%. The computers used were (i) Mac Pro with two 2.4 
GHz quad-core Intel Xeon processors and 16 GB mem- 
ory and (n) Lenovo Y570 laptop with Intel i7-2670QM 2.2 
GHz processor and 8 GB memory. With the turbo boost 
processor frequency of 3.1 GHz, the Lenovo laptop is faster 
than Mac Pro and has a higher Mathematica benchmark, 
1.1 vs 0.7. However, with the double number of processor 
cores, Mac Pro has an edge in the overall speed of our 
computations. The computation time was about two days 
for a typical relaxation curve. 




Fig. 1. Typical time dependence of spin average and the en- 
ergy bias of a magnetic molecule. Change of occurs when 
W goes through zero. 
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Fig. 2. Crossover of the spin auto-correlation function from 
weak to strong DDI that suppresses coherence of spin tunnel- 
ing. 



4 Numerical results 

The dynamics of the system in the absence of the external 
bias and nuclear spins is in accord with the tunneling- 
window concept. For any given spin, most of the time 
Wj » Z\ in Eq. ([7]), so that the spin is merely precessing 
around the z axis and no tunneling occurs. Only when 
Wi < A (dipolar bias within the tunneling window) the 
spin rotates around the x axis and CTz changes. For Eo 3> 
A tunneling events are rare and spin relaxation is very 
slow, that results in long computation times. In this limit, 
large system sizes are mandatory, because for the number 
of spins not large enough it can happen that no spin is 
within the tunneling window and the relaxation gets stuck 
forever. 

Typical time dependence of the spin polarization 
and the energy bias W oi & magnetic molecule is shown 
in Fig. [T] One can see that changes only when W ap- 
proaches zero and both time dependences are stationary 



quasi-random processes. Note that Ed is a measure of the 
dipolar field produced by one spin on its neighbors. As 
dipolar fields from different spins add up, the total bias 
W in Fig. [T] exceeds Eo by a factor of up to 10. 

Let us consider now the results for the spin autocor- 
relation function. The results for Ejj/A 1 are shown in 
Fig. [21 For a weak dipolar interaction Eo (or large tunnel 
splitting A) there are damped oscillations at frequency 
A/h showing coherence, however. Although A can be in- 
creased by a strong transverse field, Ed/ A always remains 
a large parameter. With increasing Ed/ A in Fig. [21 oscil- 
lations disappear and spin tunneling becomes incoherent. 

Figure [3] shows that for a strong dipole-dipole inter- 
action spin relaxation becomes very slow. The results of 
Fig. [3] for strong dipolar field are represented in the scaled 
form in Fig. [3] that shows a t^/^ law for the magnetization 
relaxation at short times. The results can be fitted to the 
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Fig. 3. Spin autocorrelation function for strong DDI is decay- 
ing very slowly. 




Fig. 5. Testing different exponents p for the short-time quan- 
tum spin CP, Szz{t) = 1 — at^, on a simple cubic lattice. 



(o(O)a(f)) 




Fig. 4. Spin autocorrelation function in the scaled form. Solid 
line is the stretched exponential of Eq. (|17|) . 

stretched exponential 
^..(t) = e--K.)-, ry,-j{0''- (17) 

Both the exponent in the spin relaxation function and the 
relaxation rate differ from the previously obtained result 
of Eq. ([2]). Although numerical calculations have been per- 
formed for Mni2 Ac, one can expect the same result, up 
to a numerical coefficient, also for Feg. Using parameters 
below Eq. © for Peg, one obtains r3/4 ~ 10"'* s'^. This 
by an order of magnitude exceeds the t^/^ relaxation rate 
1^1/2 ~ 10~^ s~^ measured in completely disordered Fcg 
at zero external bias. Fig. 3 of Ref. [55] . The relatively fast 
spin relaxation in the model with only dipolar interactions 
does not support the narrative of DDI blocking spin tun- 
neling so strongly that rapidly fluctuating nuclear spins 
are needed to open up a bigger tunneling window. In fact, 
this conclusion already follows from Eq. ([2]). 

As another test for the short-time dependence of the 
spin autocorrelation function, calculations for the simple 




-15 -10 -5 5 10 15 

Fig. 6. Distribution of the dipolar bias on a magnetic molecule 
in Mni2 Ac in a disordered state is nearly Gaussian with a small 
triangular distortion and extends up to ~ 20Ed- 

cubic lattice have been done. In Fig. [SJ this dependence 
is shown in three forms that allows one to choose be- 
tween three different values for the time exponent. Clearly, 
p — S/4: comes out as a winner. The same value of this 
exponent for the Mni2 Ac body-cented tetragonal lattice 
and the simple cubic lattice suggests universality of spin 
relaxation. 



5 Influence of a static external bias 

Static external bias Wq in Eq. ^ should slow down the 
relaxation because Wq has to be compensated for by the 
dipolar bias, that makes less spins being on resonance at 
any moment of time. The expected dependence of the re- 
laxation rate on Wq is given by P{Wo), the distribution 
of the dipolar bias that enters Eq. ([Ij . For a big crystal of 
Mni2 Ac in a fully disordered state P{W) shown in Fig. |6] 
is nearly Gaussian with a small triangular distortion and 
it extends up to ~ 20Ed. Here the fully disordered 
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Fig. 7. Spin autocorrelation function of Mni2 with a static 
external bias Wq. Suppression of relaxation is much stronger 
than modelled by the distribution of the dipolar bias in Fig. [S] 



state has been modelled by spin vectors cr^ (0) pointing in 
random directions, as in the dynamical calculations. 

Surprisingly, the results of numerical calculations in 
Fig. [7] show a much stronger suppression of the spin relax- 
ation by the external bias than conjectured above. Indeed, 
for the parameters in Fig. [71 Wq/A = 100 corresponds to 
Wo/{2Ed) = 5. According to Eq. ^ and Fig. the re- 
laxation rate should drop by a factor of two. However, 
the actual slow-down of spin relaxation is so strong that 
it suggests that a significant part of spins are not relax- 
ing at all, at least at the scale of these computations. For 
Wo/A = 50 (i.e., Wo/{2Ed) = 2.5) one can expect only 
a small effect of the static bias. However, the relaxation 
becomes very slow at long times. 

A likely reason for this behavior is the following. Signif- 
icant dipolar bias in Mni2 is due to regions of spins aligned 
in the same direction. Thus spins with a bias of another 
sign must be spatially well separated from the first group 
of spins. Because of this large spatial separation, tunnel- 
ing of the spins in the first group does not change much 
the dipolar bias on the spins in the second group. As the 
result, the latter have to wait a very long time until their 
dipolar bias changes via slow spatial diffusion and they 
come on resonance. 

Since the shape of the relaxation curves at nonzero 
external bias is unclear, it is difficult to extract the de- 
pendence of the relaxation rate on the bias. This problem 
requires further investigation. 

6 "Classical" spin correlation function 

As a cruder approximation, instead of Eq. ^ one could 
use the classical-like spin CF defined by 

= a,,(t)(7,^(0), (18) 

where ct^z (t) follows from Eq. ^ . Although this definition 
violates quantum properties of spins 1/2, it still captures 



Fig. 8. Testing different exponents p for the short-time "clas- 
sical" spin CF, (t) = 1 - at^, on the Mni2 lattice. 

most of the physics and is still superior to the Monte Carlo 
method as it is directly follows from the original model. 
Since the initial state is a random, the average of Sizz (0) is 
1/3, and thus one must define the average autocorrelation 
function as 

Szz{t)^^Y.^U^) (19) 

i 

that satisfies Szz{0) = 1 in the thermodynamic limit. 

As method of "classical" CFs involves only 3 dynamic 
variables per lattice site, it is two times faster than the 
method using quantum CFs. The short-time behavior of 
the "classical" CF is shown in Fig. [51 In contrast to Fig. 
[SI the results can be best fitted with p — 2/3 in Szz{t) = 
1 — at^, that differs from p = 3/4 for the quantum spin 
CF above. 



7 Influence of nuclear spins 

Magnetic atoms of some molecular magnets, such as Mni2, 
possess nuclear spins, while others such as the main iso- 
tope of Fe in Fes, do not. These nuclear spins are coupled 
to the electronic spins by the hyperfine interaction that for 
them is stronger than any other interaction. As a result, 
nuclear spins are precessing in the effective field created 
by the electronic spins and produce a random but static 
bias on the latter. This can essentially modify the Landau- 
Zener effect [44145] but does not provide a fiuctuating bias 
that could extend the tunneling window. 

On the other hand, all magnetic molecules contain a 
large number of protons (in the hydrogen atoms of their 
ligands) whose nuclear spins interact with electronic spins 
by the nuclear dipole-dipole interaction. As in the case of 
the DDI, the only realistic way to consider the system's 
dynamics is to make decoupling of quantum correlators 
similar to Eq. ([S]). After that nuclear spins can be consid- 
ered as time-dependent bias on electronic spins. 

Since protons have nuclear spin / — 1/2, there are no 
quadrupolar terms in their Hamiltonian, and their dynam- 
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ics is precession in the magnetic field B, 



i ^ ifip/h) [I X B] 



(20) 



with the frequency ujp — {fip/h)B. Here /ip is the magnetic 
moment of the proton and /i„ is the nuclear magneton. 
If no external field is applied, the dominating magnetic 
field is dipolar field from the electronic spins. While the 
fully ordered state of Mni2 Ac creates the dipolar field 
of 53 mT, disordered states of Mni2 Ac and Fes create 
a somewhat smaller field. With B = 16 mT one obtains 
ojp = 1.3 X 10^ s~^ that largely exceeds the tunneling 
frequency A/h. For Feg in zero field A/h = 1.3 x , 
so that hu!p/A ~ 10^. If a strong transverse field B± is 
applied to increase A, precession of nuclear spins becomes 
regular with the fixed frequency Up = {iip/h)B±. For Fes, 
the ratio hup/A initially strongly increases with B± but 
then decreases as A begins to increase as a high power of 
B± [33] (see experimental data in Ref. [12 )• In all cases 
hup/A is a large parameter. 

To estimate the fiuctuating bias from the protons, one 
can use the data for the tunneling resonance spread of 
approximately Bp ~ 0.5 mT (data for Fes from Ref. [26j). 
As follows from the properties of the dipolar interaction, a 
part of this bias is static and thus irrelevant here, whereas 
its another part is precessing with a frequency ujp. The 
latter creates a time-dependent bias 

W^^''\t) = Wp cos {ujpt + if,), Wp = 2Sg^iBBp. (21) 

In contrast to the model used in Ref. [22], relevant dy- 
namics of nuclear spins is precession, rather than their 
dephasing at rate T2~^. Here, one could take into account 
distribution of Wp values, as well as possible distribution 
of ujp. However, this does not bring any qualitative change 
of the results because the dominating DDI already creates 
enough randomness for a complete decoherence. For this 
reason and the sake of simplicity, numerical calculations 
have been performed for the fixed ratio Wp/ {2Eu) = 0.1 
corresponding to molecular magnets and different fixed 
values of the ratio hujp/A. Note that the nuclear bias Wp 
is by a factor of order 10^ smaller than the dipolar bias, 
the distribution of which is shown in Fig. [51 

Let us check whether the bias created by precessing 
nuclear spins of protons is fast, as required [22], by con- 
sidering Landau-Zener transitions of electronic spins. The 
typical value of the LZ parameter 



ttA^ 
2hv 



(22) 



where v is the energy sweep rate. With v = dwj;^\t)/dt 



Wp(jjp it can be estimated as 



2fujjpWp 



(23) 



Since both Wp/A » 1 and hwp/A > 1, one has e ^ 1 
and the sweep is indeed fast. 

Numerically, it is difficult to perform calculations for 
a very rapidly oscillating nuclear bias, Eq. (j2ip . with the 
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Fig. 9. Spin correlation function taking into account nuclear 
spins. The latter speed up the relaxation by a factor of about 
3 in the limit of rapid precession. 



same method as above because it requires integration of 
equations of motion with a very small time step over a long 
time. Whereas an appropriate semi-analytical method will 
be considered elsewhere, here we just perform numerical 
calculations for moderately high values of fnUp/A to see 
the effect of nuclear spins. It can be expected that in the 
case of a slow nuclear precession, £ ^ 1, the effect of nu- 
clear spins is not strong. One cannot speak of widening 
of the tunneling window because some spins are getting 
a chance to tunnel because of Wp while other spins lose 
their chance. Thus, the tunneling chance is merely redis- 
tributed between spins but no additional tunneling chance 
is created. 

On the other hand, one can speak of forced spin transi- 
tions via the Landau-Zener (LZ) effect. It has been shown 
that interactions, including the DDI, severely suppress LZ 
transitions [S]. For a fast sweep, however, interactions do 
not play a big role, and the standard LZ effect with a small 
transition probability takes place. In spite of a small tran- 
sition probability at fast sweep, there are repeated level 
crossings for an oscillating bias, and the effect of many 
crossings accumulate. It was shown that a combination 
of a time-linear sweep and a rapidly oscillating sweep re- 
sults in the same transition probability as the pure linear 
sweep [33 • Thus one can expect a more substantial effect 
of nuclear spins on low-temperature relaxation in molecu- 
lar magnets in the realistic case of rapidly precessing nu- 
clear spins. This effect should saturate at high precession 
rates of nuclear spins. 

The results of numerical calculations in Fig. [H] confirm 
the theoretical conjectures made above. In the case of a 
slow precession hujp/A — 0.1 (thus e ~ 1) the effect of 
nuclear spins is small. For a fast precession hiOp/A ~ 10 
(thus £ ~ 0.01) it is substantial, creating a speed-up of 
relaxation by a factor of about 3. However, such a speed- 
up is not a profound effect that could be expected from the 
argument that tunneling window A is replaced by that of 
Wp ^ Z\. In fact, this effect should be somewhat smaller 
because only a part of the nuclear bias is oscillating while 
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another part is static. On the other hand, spin relaxation 
in the absence of nuclear spins, Eq. (IT7| . is already fast 
enough, so that one really does not need an additional 
speed-up from nuclear spins to describe experiments. 

It is difficult to further increase within the present 
computational scheme and time limitations. However, e ~ 
0.01 for hujp/A = 10 means already a pretty fast sweep, 
and the results can be expected to saturate in the fast- 
sweep limit. Indeed, calculation for hujp/A — 100 with 
a 10 times smaller integration step and 10 times shorter 
time interval (to preserve the computing time) shows the 
same results as for hwp/A = 10. 



8 Discussion 

Results for the dipolar-controlled low-temperature spin 
relaxation in molecular magnets in the disordered state 
presented above have been obtained with a direct method 
based on the actual quantum spin dynamics of the system. 
An essential approximation made is decoupling of quan- 
tum correlations at different sites, the quantum mean field 
or Hartree approximation. Although the results of this 
approximation may differ from the exact solution of the 
Schrodinger equation (in particular, for the Landau-Zener 
effect in systems with DDI [H]), one can expect that this 
approximation captures the essential physics and its re- 
sults are qualitatively correct. In any case, this approach 
is more fundamental than previous Monte Carlo simula- 
tions and yields the relaxation curves in terms of the real 
time rather than of Monte Carlo steps. 

Most of the results have been obtained for Mni2 Ac lat- 
tice in the case of a pure DDI without nuclear spins. For 
other molecular magnets such as Feg the results should be 
qualitatively the same. The results show that relaxation 
is due to rare tunneling events in the actual case of ran- 
dom dipolar bias much greater than tunnel splitting. The 
relaxation law is a stretched exponential in the range of 
not too long times. The exponent in the stretched expo- 
nential extracted from fits to the data obtained is p = 3/4 
for both Mni2 Ac and simple cubic lattices, differing from 
the experimental law t^/"^ |25l26j and Monte Carlo results 
of Refs. [29130131] . On the other hand, it is in a quali- 
tative accord with Monte Carlo results of Refs. [32 33], 
p = 0.7 — 0.73 for a face-centered lattice. 

The effect of nuclear spins has been included on the 
basis of a simplified microscopic model and it has been 
shown that nuclear spins speed-up the relaxation in the 
realistic case of their fast precession. However, the ob- 
served speed-up is moderate and it does not support the 
idea of nuclear spins opening a huge tunneling window. 

The deviation of the relaxation law obtained above 
from the experimentally observed may be a consequence 
of either thermally assisted tunneling or correlations ac- 
companying dipolar ordering, or both. The problem is that 
to get rid of the ordering effects, the temperature has to 
be above 1 K. However, at such temperatures thermally- 
assisted tunneling should be already strong. Although pop- 
ulations of all states except of the ground-state doublet 



are still negligible, there is a competition between differ- 
ent ways to cross the barrier, pure tunneling considered 
here and thermally-assisted tunneling over higher dou- 
blets. Since the tunnel splitting strongly increases with 
the energy, thermally- assisted tunneling remains competi- 
tive down to the temperatures below 1 K |49I50| . Thus the 
theory developed above may have no applicability range, 
strictly speaking. Still it makes sense as the most basic 
microscopic theory of dipolar controlled spin tunneling. 
It should be noted that none of the preceding theoreti- 
cal works on low-temperature spin relaxation dealt with 
these two effects, dipolar ordering and thermally assisted 
tunneling. The only investigation of the interplay of the 
dipolar ordering and dipolar controlled relaxation is that 
in the recent experimental work on Er, Ref. [510. 

The discussion above shows the tasks of future inves- 
tigations. It would be interesting to extend the method to 
higher temperatures into the regime of thermally assisted 
tunneling, as well as to lower temperatures into the re- 
gion of dipolar ordering. The latter is challenging because 
it requires preparation of the initial state having a partic- 
ular dipolar energy corresponding to a given temperature. 
As at finite temperatures there are spin-spin correlations, 
one cannot use the autocorrelations function. As a result, 
self-averaging does not occur and one has to make averag- 
ing over many initial states, even for a system of a large 
number of spins. 
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Appendix: Spin relaxation function via spin cor- 
relation function 

Let us define the linear response of the normalized pseu- 
dospin with respect to the energy variable ^ = Sg^sBz- 
The linear response has the form 

/oo 
dt'L{t-t')at'). (24) 
-oo 

where, because of the causality, L{t) ~ for r < 0. We 
are interested in the response to the step ^(t) — fo^(~i) 
(field switched off at t = 0) that reads 

(a.),=?o / dt'L{t-t'), t>0. (25) 

It is convenient to express this response in terms of the 
dynamic susceptibility. Response to the oscillating field 

{{t) = Coe"'"^* has the form 

{<^z)t - eo / dt'L{t - t')e-'^'' = ^oe-*"*x(^), (26) 
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where 



dTL{T)e'' 



(27) 



is the susceptibility with respect to ^. Inverting this for- 
mula, one obtains 



L{r) ^ r pe-— 
Now Eq. ([^5]) becomes 

/•OO 



X(w), 



-0. 



(28) 



dui 



2iT iuj + e 



Xi^)- (29) 



The imaginary part of the susceptibility is related to 
the spin-spin CF by the fluctuation-dissipation relation 



dtcos{ujt)S^^{t). (30) 



To calculate the step response, Eq. ([^ , one needs the full 
susceptibility, and x'(w) can be found from the Kramers- 
Kronig relation 



1 f y" 
X'(^) = -r / du:'^ 

TT J-oo ^ 



Substituting Eq. ([50j) into Eq. (PT|) one obtains 

/OO 
-OO 

where 



(31) 



(32) 



duj' , hu' 
e^"*tanh— - 

LU' — U! Kb! 

^e'("+-)*tanh^(" + "\ 



(33) 



This can be rewritten in the form with no singularity in 
the integrand at z — > 0: 



.».t.„„u ^(w + z) 



e*^' tanh ■ 



--tanh^("-") 



(34) 



At high temperatures one can expand the tanh terms and 
after integration obtain 



1 



e*"* [i7ra;sign(t) + 27r(5(t)] . (35) 



Inserting this into Eq. (15^ and taking into account that 
Szz(t) is an even function, one obtains 



1 — uj dt sm{ujt)Szzit) 



Combining this with Eq. pop at high temperatures 
2io 



X"{^) = 



knT 



B-L Jq 



dtcos{u}t)Szzit) 



(36) 



(37) 



yields 



x(^) 



knT 



l + iu I dte^^'Szzit) 





(38) 



Now for the step response from Eq. follows 
2^0 



ksT 7-00 27r icij + e 



l + iu dt'e"^''Szz{t') 
Jo 

(39) 

Here the first term vanishes and in the second term one 
can set e = 0. This yields 



/cq 



Szz{t), 



(40) 



where {crz)cq = 2^0/(^5^) is the equilibrium spin polar- 



ization. 
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